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ABSTRACT 


The Space Shuttle program has relied heavily on simulation 
throughout all phases of development and operation. Real-time, 
man-i n-the-1 oop simulation has served the NASA manned-space-flight 
program by providing the means to evaluate systems design and integrated 
systems performance in a simulated flight environment as well as provide 
a means to train flight crews. 

New challenges are presented by the development and operation of a 
permanently manned Space Station. The assembly of the Space Station, 
the transferral of payloads and the use of the Space Station manipulator 
to berth the Orbiter are operations I critical to the success of the Space 
Station. All these operati ons |are examples of constrained motion among 
the bodies associated with the Orbiter and Space Station system. 

Current real-time simulations of the Orbiter and Space Station do 
not have the capability to model the dynamics of constrained motion. 
Determining an efficient, yet general, method for constrained motion in 
multibody systems is one essential key to high-fidelity of these 
critical operations. 

This report described the state-of-the-art of formulating the 
governing dynamical equations of motion for constrained systems. The 
uses of the two basic problems in multibody dynamics are discussed. The 
most efficient formulations of the equations of motion are addressed 
from the point of view of completeness. The issues surrounding 
incorporating the constraints into the equations of motion are 
presented. Finally, an overview of the strengths and weaknesses of the 
current methods is given, along with some recommendations for further 
research. 


INTRODUCTION 


The planned scenarios involving the Shuttle Orbiter, the Space 
Station and their respective remote manipulators are examples of complex 
dynamical systems subjected to kinematic constraints, see Figure 1. 



i 

Figure 1 - Constrained motion between TRS and satellite 

In order to develop realistic models of these activities, some of 
the bodies may be considered to be rigid and others flexible. An 
analysis of the flexibility of the bodies requires a complete knowledge 
of the system dynamics so that quantitative expressions for the joint 
forces and moments may be obtained. Having the capacity to formulate 
completely the equations of motion with the appropriate constraint 

equations and solve them is essential to any simulation of these ! 

systems. 

NASA has relied heavily on simulation throughout the Space Shuttle j 

program's development and operation. Real-time, man-i n-the-1 oop 

simulation has provided the means to evaluate systems design and ! 

integrated systems performance in a simulated flight environment and to 
train flight crews. 
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Current real-time simulations of the Orbiter and Space Station do 
not have the capability to model the dynamics of constrained motion. 
Determining an efficient, yet general, method for simulating constrained 
motion in multibody systems would enhance the fidelity between the 
man-in-the-loop simulation and critical on-orbit operations. 

This report addresses the efficiency and completeness of current 
formulations of the dynamics of constrained systems. Efficient 
formulations are necessary since the goal is to run the simulation in 
real-time. The formulation should also be complete since flexible 
bodies may be included. 

The balance of this report is divided into three parts with the 
following part providing a discussion of the equations of motion. The 
next part discusses incorporating the constraint equations. Concluding 
remarks are given in the last part. 


EQUATIONS OF MOTION 


There are two basic problems in multibody dynamics. The first is 
the problem of determining the motion of a system from a set of applied 
forces; it is referred to as the forward dynamics problem. The second 
is the problem of determining a set of forces required to produce a 
prescribed motion in a system; it is referred to as the inverse dynamics 
problem. 

A solution of the inverse dynamics is essential in the dynamic 
control of systems and many efficient formulations have been developed 
to implement real-time dynamic control. A solution of the forward 
dynamics is essential in the simulation of dynamic systems. Usually the 
forward dynamics are found by first formulating the inverse dynamics and 
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then solving for the motion. Recently, several methods for improving 
the efficiency of the forward dynamics solution have been developed. 


The three main approaches towards deriving the dynamic equations of 
motion for multibody systems with rigid elements have been the 
Newton-Euler, the Lagrange and the Kane methods. In all these methods 
efficiency is derived from the structure of the computation. P ir$t the 
linear and angular velocities and accelerations are computed recursively 
from the reference body to the end of the chain, then the forces and 
torques are computed recursively from the end to the reference body. 



Figure 2 - A general chain system 


Consider a multibody system such as shown in Figure 2. An 
accounting system for the system connectivity may be developed by 
arbitrarily selecting one of the bodies as a reference body and calling 
it Next, number the other bodies of the system arbitrarily. This 

numbering of the bodies can be used to describe the chain structure or 
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topology through the "body connection array" as follows: Let L(k), k = 
1, ...,N be the array of the adjoining lower body of B^. The lower body 
is defined as the body adjoining B k that is closest to the reference 
body [1]. 


The angular velocity of a typical body B k in the inertial frame, 
B q , is readily obtained by the recursion formula 


k 

U) 


. L(k) ♦ L(k) c/ 


( 1 ) 


J k and B L(k) 


where M k and W L ^ are the angular velocities of B, 

^ 1 f k ^ k * 

to the inertial frame and v is the angular velocity of B^ 

to B L(k) , 


relative 
relati ve 


The angular acceleration of B k in the inertial frame may be 
obtained by differentiating equation (1) in the inertial frame, leading 
to the recursion formula 


L(k) , L(k)^k L(k) x L(k) u k (2) 

— k — — — — 

where °c k and « L ^ k ^ are the angular accelerations of B K and B. , . ) 

— — i ii,\ I, 

relative to the inertial frame and v is the angular acceleration of 
B k relative to B^ k ). 

Recursion formula can also be developed for the position, velocity 
and acceleration of a point on any link [1-4]. 


A chain system of N rigid bodies will, in general, have 6N or fewer 
degrees of freedom. This comes as a result of the fact that each rigid 
body can have up to 6 degrees of freedom. Kinematic constraints between 
members of the system reduce the number of degrees of freedom through 
contact between kinematic pairs. Usually, the analysis is simplified by 
introducing the kinematic constraints due to joints connecting the 
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bodies of the system at an early stage. The systems' s configuration is 
then said to be completely determined by n generalized coordinates, 
where n is the number of degrees of freedom of the system. However, it 
may occur that a contact between bodies may be made or broken as the 
motion of the system progresses. Figure 3 shows a case where two more 
bodies of the system in Figure 2 are in contact. In order to 
distinguish between the two types of constraints the former will be 
called a permanent constraint and the latter will be called a temporary 
constraint [5], 



Figure 3 - Temporary constraint between bodies 4 and 6 


Many procedures have been, developed for both the inverse and 
forward dynamics of systems with permanent constraints using recursive 
relationships for efficiency. References [2, 3, 6] represent recently 
reported research efforts in rigid body inverse dynamics using recursive 
formulations of the Newton-Euler , Kane and Lagrange methods. 
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respectively. Recursive formulations of the forward dynamics of rigid 
body systems with permanent constraints are significantly -fewer in 
number and more intricate in their derivation. References [7-9] 
represent the currently reported procedures of this type. 

Armstrong [7] hypothesized the existence of a linear recursive 
relationship between the motion of a forces applied to a body and the 
motion of and forces applied to its neighbors. He defines a set of 
recursion coefficients for each body and shows how the coefficients for 
a body may be calculated in terms of those of one of its neighbors. - The 
coefficients are then used to calculate the acceleration of each body. 
The computational complexity of this method is 0(n), but the coefficient 
of n is quite large. The method is applicable to robots with spherical 
joints, but a modification is included in Appendix II [7] for revolute 
joints. This modification increases the computational requirement 
significantly, although the methods remains 0(n). 

Walker and Orin [8] describe four methods of solving the forward 
dynamics problem for systems with revolute or prismatic joints. Three 
of these methods are based efficient techniques for solving the inverse 
dynamics problem. The computational complexity of each of these methods 

O 

is 0(n ), however the coefficients of the higher powers of n are 

relatively small. The fourth method discussed is based upon a recursive 

technique for constructing the moment of inertia matrix. The concept of 

composite center of mass and moment of inertia matrix of a partial set 

of bodies at the end of the chain is used in this recursion. The fourth 

method has a computational complexity of 0(n ), however the coefficient 
2 

of the n term is large. 

3 

The basic approach of Walker and Orin's 0(n ) algorithm can be 
summarized by following Featherstone [9], The equations of motion can 
be expressed in the form 
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(3) 


T = A(C[) £ + C(£, £, X) 

where t_ is the vector of generalized forces acting at the joints; 

£, £ are the generalized coordinate vector and its time derivatives; 

A is the generalized inertia matrix and c is the vector of Coriolis, 
centrifugal and external forces x.. 

An algorithm for solving inverse dynamics can be thought of as an 
implementation of a vector function in the form 

t = d(£, £, £, x) (4) 

which states _t is obtained from the generalized coordinates, the 
generalized speeds, the generalized acceleration and the external force 
acting on the system. It is seen that (cf (3)) 

cCg.. £, x) = d(£, £, 0, x) (5) 

the calculation of £ is the special case of calculating the inverse 
dynamics when the generalized acceleration is zero. 

The forward dynamics calculation is conveniently partitioned into 
three steps: the calculation of c_, the calculation of A, and the 

solution of the set of simultaneous equations 

Ac[ = _r - £ (6 ) 

For £. The computational complexity of the first step is 0(n), the 
2 3 

second step is 0(n ) and that of the third is 0 ( n ). This approach is 
referred to as the composite rigid-body method. 

Featherstone [9] describes an extension of Armstrong's method to 
system with revolute and prismatic joints and uses a spatial notation 
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consisting of 6X1 spatial vectors, 6X6 spatial transformation matrices 
and 6X6 spatial inertias. He refers to this method as the articulated- 
body method. 

The articulated-body method starts with the equation 

q = _h(<l, l - c) (7) 

which states that | is a function of £ and x - £ without implying the 
need to calculate a generalized inertia matrix. 

The calculation is again conveniently partitioned into three steps: 
the calculation of c, the calculation of the inhomogeneous articulated- 
body inertias for each body starting with the end effector and 
proceeding sequentially to the reference body, and the generalized 
accelerations are then calculated for each body starting at the 
reference body and proceeding sequentially to the end effector. The 
computational complexity of each step is 0(n). 

The computational complexity of Armstrong's method is greater than 

3 

either of the most efficient of Walker and Orin's 0(n ) methods or 
Featherstone's articulated-body method. The best Walker and Orin method 
is more efficient than the Featherstone for values of n less than twelve 
[9], Either the Armstrong nor the Walker and Orin nor the Featherstone 
methods provide for the calculation of the constraint forces of the 
joint. In this sense these algorithms are not complete. The next part 
describes some methods of incorporating the constraint equations and 
finding the'constraint forces. 


INCORPORATING CONSTRAINTS 

References [5, 10-15] represent a partial list of recently reported 
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research efforts on dynamics formulations for constrained motion. The 
equations of motion for a system of rigid bodies with open chain or tree 
structure (as in Figure 2) may be written as in equation (6). In 
systems where additional restrictions have been placed on the motion of 
the bodies, these constraints can usually be placed in the form 

Bi = g (8) 

where B is an mxn matrix with m<n and £ is n vector. Both B and £ are 
functions of £ and time t. Holonomic constraints can always be written 
in the form of equation (8). Simple or Pfaffian nonhoinonomic 
constraints can also be written in the form of equation (8). 

There are two fundamental approaches to solving the equations (6) 
and (8) for such systems [14]. One can introduce unknown generalized 
forces of constraint between the constrained bodies by the technique of 
undetermined multipliers and then solve the equations of motion 
simultaneously with the constraint equations to determine the 
generalized constraint forces as well as the kinematic variables. 
Alternatively, one can use the constraint equations to reduce the 
dimension of the equations of motion to be solved, eliminating the need 
to represent or solve for generalized constraint forces. The second 
approach has the advantage of providing, for solution, a set of 
equations of minimum dimension. The difficulty lies in its 
impl ementation. 

In the special case of holonomic constraints, one might integrate 
equation (8) and solve explicitly for a subset m of the generalized 
coordinates in terms of the remaining 1 = n - m. Substitution for these 
generalized coordinates in equation (6) permits reduction of the set of 
equations to be solved to 1 second order equations. In general, this 
reduction is not easily accomplished. 
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If equation (8) is nonintegrable (i.e., nonhol onomic) , then in of 
the generalized coordinates cannot be solved for in terms of the 
remaining 1. One would like to solve for m of the elements of £, but 
the result cannot be substituted effectively into the second order 
equation (6). This obstacle can be avoided by converting equation (6) 
to first order form [10]. A large obstacle is the potential for the 
rank r of the matrix B in equation (8) to be less than the row dimension 
m. These state of affairs could occur in a system changing 
configuration, resulting in singularities in the numerical solution. 
Other methods which avoid all singularity problems exist and are 
described next. 

The n generalized speeds, £, may be expressed in terms of 1 new 
generalized speeds, £, by hypothesizing a relationship of the form 

£ = &u + 6 (9) 

where 8 is an nxl matrix and { is an n vector. 

The solution to equation (8) for £ falls into one of two 

categories. If £ does not lie in the column space of B (i.e., if £ 

cannot be formed by a linear combination of the columns of B) , the 
constraints are inconsistent, the problem is ill-posed and no solution 
can be found for £. If £ does lie in the column space of B, at least 
one solution exists for £. Let Z be an nxl matrix with linearly 

independent columns that form a basis for the null space of B, that is, 

BZ = 0 (10) 


and let Y be an nxm matrix whose columns complete a basis for the vector 
space R n . The general soluation to equation (8) is 

£ = Z£ - Yy (11) 
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where is a vector of m unique scalars such that 


I - BY i (12) 

and 2 is a vector of 1 independent quantities describing the constrained 
system. The vector u^ of generalized speeds for the constrained system 
may be defined as 

u_ - £ (13) 

By substituting from equation (13) into (11), one obtains an 
equation of the form of equation (9) with 


$= Z 


(14) 


and 


6 = -Yy (15) 

The preceding development demonstrates that it is always possible 
to find an equation of the form of equation (9) for a system subjected 
to constraints in the form of equation (8). However, it does not give 
an operational procedure for producing equations (14) and (15). The 
matrices Z and Y are not unique, and many approaches for constructing 
them can be envisioned. The problem is equivalent to that of 
constructing a generalized inverse of B, of which the Moore-Penrose 
generalized inverse is a particular example. Singular value 
decomposition [12, 14] and the zero eigenvalue theorem [11], are two 
methods of calculating this inverse. The remainder of this part will 
briefly outline two methods using Kane's equations of motion. 

Kane's equations of motion [10] may be written in the remarkably 
simple form 
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f + f* = 0 


( 16 ) 


where _f is the n-vector of generalized active forces and f_* is the n 
vector of generalized inertia forces. The set of equations (6) may be 
efficiently cerived from the set of equations (16) by relatively simple 
scalar products and recursive kinematic relationships [3, 4], 

Kane's equations for constrained systems may be found by solving 
equation (8) for m of the £ in terms of the remain 1, as described 
before. This leads to a reduced set of equations of motion. 

k + k * = 0 (17) 

where k and k* are 1 -vectors of reduced generalized active and inertia 
forces. Equations (17) together with equations (3) then constitute a 
system of n equations for the n-vector £. 

For complex multibody systems the solution of equations (3) for m 
of the £ may not be convenient or even possible. Also, the generalized 
constraint forces have been eliminated from the analysis. Although this 
reduces the dimension of the problem, it also means the analysis is 
incompl ete. 

To develop a method which includes the generalized constraint 

equations and retains the advantageous features of the previous method, 

let the constraining forces and moments associated with the constraint 

equations (8) be represented by a series of forces f'. applied at points 

P., i = 1, . . . , N ' , together with a series of couples with moments m 1 . 

applied to bodies B., j = 1, ...,M'. Assuming that m < 3 (N‘ + M'), the 

J 

following analysis is applicable [15], 

Let the P. and B. have specified motions given by Vp.(t) and u>. ( t). 
Let the scalar components of these vectors, relative to a set of basis 
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vector fixed in the reference body, be collected into arrays and w 
whose elements may be identified with the elements of £ in equation (8). 
That is. 


O) 

i 


9i 


m 

= E 


j=l 


b . . q . 
ij H J 


(18) 


for i = 1,...,3N' for the V. and i = 3N‘ + 1, ...,3(N'+M‘) for the w. . 


Next, let the components of £' . and M' . in the directions of the 

* J 

components of Vp- and w. be designated as 4> . and y.. Then the power, P, 
produced by the constraining forces and moments is 

3N' 3M‘ 

P = E 4). V. E y uk (19) 

1=1 j=l 

Let A be an m-vector whose elements are <J>^ and y . P may now be written 
as 

P = A T £ = ^ TB i (20) 

where the superscript T designates the transpose. 


Since neither _f ' . nor mV are functions of £, the generalized 
forces associated with these constraining forces and moments are 

1' = 9P/8q = B T X (21) 

where _f' is the generalized force vector. 

In view of the foregoing analysis, let the generalized active 
forces be separated into two parts consisting of: those developed from 

constraint forces and moments and those developed from applied forces 
and moments. Let these be designated _f' and f_. Then Kane's equations 
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are 


_f ' + f_ + 1* = 0 (22) 

Substitute for f 1 in equation (22) using equation (21) 

rT A + 1 + 1* = 0 (23) 

Taken together with the constraint equations (8) they constitute a set 
of n+m equations for the n elements of £ and the m elements of K 

Suppose that an nxl matrix has been found such its columns form a 
basis for the null space of B. Then by premultiplying equation (23) by 
we obtain 

e T f + B T f* = 0 (24) 


si nee 


86 = B T B T = 0 (25) 

Equation (24) together with equations (8) then represent a set of n 
equations for the n elanents of £. Once £ and £ are determined, A may 
be determined by back substitution into equation (23). 

The method may be viewed as a generalization of the psuedo- 
coordinate method of Kane [10], The analysis is developed with the 
system unrestrained. However, a difference is that the constraint 
forces and moments are formally introduced through a product of 
undetermined multipliers with coefficients of the constraint equations. 

The concept is simple, but the implications extend beyond the task 
of finding constraint forces and moment components. For example, when 
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compared with Lagrange's equations with undetermined multipliers, the 
procedures are similar. However, with Kane's equations the generalized 
speeds may be employed as the fundamental kinematic variables. Also, 
holonomic as well as simple nonholonomic constraints may be accommodated 
directly. Finally, the constraint force and moment components may be 
either evaluated or eliminated from the analysis. In the latter case 
the resulting equations are equivalent to the standard form of Kane's 
equations for constrained systems. 

To illustrate these concepts consider Figure 4 which depicts the 
generalized forces in the n-dimensional space of the generalized speeds. 
The force triangle is a representation of equation (22). In this 
context, the generalized forces k and k*, used in Kane's equations for 
constrained systems, are seen to be projections of _f and f* onto the 
directions the generalized inverse of B. 



Figure 4 - A geometric interpretation of the generalized >rces 
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CONCLUSIONS 


Very efficient procedures have been developed for solving the 
inverse dynamics problem. Representative solution method using 
recursive formulatins of the Newton-Euler, Kane and Lagrange methods are 
described. A generalization of Kane's method [15] has been shown to be 
capable of efficiently and completely accommodating various constraints 
applied to a system whereas the others have not. This is not to say 
that this cannot be done. More research is certainly appropriate to 
develop parallel techniques using the Newton-Euler and Lagrange methods. 

The forward dynamics problem has been more elusive. Although very 
efficient procedures have been developed for the special cases of 
permanent constraints of the revolute, prismatic and spherical types, a 
corresponding generalization has yet to be developed. Indeed, none of 
the forward dynamic algorithms available can efficiently computer the 
constraint forces and moments for even permanent constraints. Research 
is to determine efficient and complete methods of generalizing the 
solution to the forward dynamics problem, if this is possible. 
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